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Abstract. The equations governing null and timelike geodesies are derived within 
the 3+1 formalism of general relativity. In addition to the particle's position, they 
encompass an evolution equation for the particle's energy leading to a 3+1 expression 
of the redshift factor for photons. An important application is the computation of 
images and spectra in spacetimes arising from numerical relativity, via the ray-tracing 
technique. This is illustrated here by images of numerically computed stationary 
neutron stars and dynamical neutron stars collapsing to a black hole. 

PACS numbers: 04.25.D-, 95.30.Sf 
1. Introduction 

The computation of trajectories of photons or test-mass particles in the Kerr metric 
is a topic of major importance in relativistic astrophysics. This notably allows the 
investigation of spacetime properties around black holes (see e.g. [U El El IH El E] and 
references therein), the aim being to determine the black hole's mass and spin and to 
test general relativity (GR). Photons and test-mass particles follow spacetime null and 
timelike geodesies, respectively. Their motion is thus governed by the so-called geodesic 
equation. 

However, within the framework of metric theories, strong tests of GR require 
to compare the Kerr geometry with geometries generated by alternative models of 
compact objects. The metric is then generally not known analytically and must be 
computed numerically. Rotating gravastars and boson stars are examples of such 
objects. Numerical metrics being almost exclusively computed using the 3+1 formalism 
of GR (see e.g. |7]), it is quite useful to derive the geodesic equation within this 
framework. This is a way to build an optimized ray-tracing algorithm. 

In addition to the GR tests around astrophysical black holes, another field 
of application is the visualization of computer-generated spacetimes, resulting from 
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numerical relativity studies of sources of gravitational radiation, such as gravitational 
collapse or coalescing binary compact objects [8j [9]. Such spacetimes are generally 
computed with the 3+1 formalism, and this motivates the design of a ray-tracing 
algorithm based on a 3+1 geodesic equation. 

A ray-tracing code capable of using a 3+1 metric, GYOTO, has recently been 
developed in our group [10]. This code is written in C++, is open source and can 
be freely downloaded from [TJJ. It computes null and timelike geodesies, both in the 
Kerr metric and in any numerically computed spacetime. The goal of this article is 
to derive the 3+1 geodesic equation that allows GYOTO to compute images and spectra 
in numerically generated spacetimes, and to give the first examples of astrophysical 
interest of this capability. To our knowledge, in previous works, the geodesic equation 
has only been integrated in numerical spacetimes for the purpose of locating event 
horizons [121 E2 EL but not to form images nor to compute spectra. 

The plan of the article is as follows. Section [2] derives the 3+1 geodesic equation. 
Section [3] derives the 3+1 expression of the redshift factor, useful for ray-tracing 
computations. Section [4] presents the first applications of ray-tracing in numerical 
spacetimes considering stationary and collapsing neutron star spacetimes. Finally, 
Section [5] gives conclusions and perspectives for future works. 

2. 3+1 geodesic equation 

2.1. Framework 

Let (Ai,g a p) be a 4-dimensional spacetime, i.e. a 4-dimensional smooth manifold M. 
endowed with a pseudo-Riemannian metric g a p, of signature (—,+,+,+). We denote 
by V a the Levi-Civita connection associated with g a p. The 3+1 formalism of GR (see 
e-g- pU El EE]) is based on the assumption that (M., g a p) is globally hyperbolic, so that 
it can be foliated by a one-parameter family of spacelike hypersurfaces (E t ) teR . Let 
n a be the future-directed unit normal to the hypersurface Ej. n a is collinear to the 
gradient of t, the proportionality factor defining the lapse function N: n a = —NV a t. 
The unit timelike vector n a is the 4-velocity of the so-called Eulerian observers Oe, i.e. 
the observers whose worldlines are orthogonal to the hypersurfaces E t . 

Using standard notations, we denote by 7 Q( g the metric induced by g a p on each 
hypersurface E t (first fundamental form). Since Et is assumed to be spacelike, 7 Qj g is a 
Riemannian metric (i.e. positive definite). One has 



and 7°^ is the orthogonal projector onto Ej. We denote by D a the Levi-Civita connection 
associated with the metric 7 a/ 3 on E t . The 4-acceleration of an Eulerian observer is 
a a := n M V^n a and obeys 



lap = 9ap + n a n/3 



(1) 



a Q = D a In N. 



(2) 



In particular, n^a 11 = 0, so that a a is tangent to E t . 
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The extrinsic curvature tensor, or second fundamental form, of the hypersurface E t 
is defined by 

Kap ■= ~Y a r^,n u . (3) 

One has K^n 11 = as well as the useful relation (see e.g. [7j) 

V/3n a = -K af3 - D a In N n p . (4) 

In this article, we shall consider only coordinate systems on M. that are adapted 
to the 3+1 foliation (£ t ) feK , i.e. coordinate system (x a ) such that x° = t. The three 
remaining coordinate^] (x l ) span the hypersurfaces S t : by construction the vectors 
d/dx l are tangent to S t . The vector d/dt is transverse to S 4 and its 3+1 decomposition 
defines the shift vector f3 a : 

-j = Nn a + f3 a , with V^ = 0. (5) 

The knowledge of the lapse function N , the shift vectoij|] f3 % and the induced metric 7^ 
is sufficient to reconstruct the spacetime metric g a p according to 

g^ dx» dx u = -N 2 dt 2 + 7ij (dx* + /?Mt) (dx j + /3 j dt) . (6) 

In the 3+1 formalism, the solution of the Einstein equation is obtained by solving a 
Cauchy problem (with constraints) for (7^-, Ify): some initial values being given on a 
hypersurface S , obeying the constraint equations, the fundamental forms (•y i j,K i j) are 
evolved in time to construct the whole spacetime [Zl El E] - In this formulation, the lapse 
and shift are not dynamical variables; their role is to set the coordinates. 

2.2. Geodesic equation in 3+1 covariant form 

Let us consider a particle V of 4- momentum p a . V can either be a photon, in which 
case p^ = 0, or a massive particle, of mass m = yj—p^p* 1 . If V is subject only to 
the gravitational field, its worldline £ is a either a null (photon), or a timelike (massive 
particle) geodesic of (M.,g a p). The 4-momentum then obeys 

P^^p a = 0. (7) 

This is the geodesic equation in covariant 4-dimensional form. To write it in 3+1 form, 
we start by performing an orthogonal decomposition of p a according to 

p a = E(n a + V a ), with n/" = 0. (8) 

The scalar E is nothing but the energy of V as measured by the Eulerian observer O^', 
indeed, n a is the 4- velocity of Oe, and then (J8J) implies E = —p^n* 1 . The vector V a is 
by construction tangent to S 4 and coincides with the 3-velocity of V as measured by 

% Latin indices span {1, 2, 3}, whereas Greek indices span {0, 1, 2, 3}. 

§ We write j3 l when we consider the shift as a tangent vector on the manifold E t and /3 a when we 
consider it as a vector onM, as in ([5]); then (3° = 0. We extend this notation to all tensor fields tangent 
to E t , in the sense that their contraction with n a or n a vanishes (for instance "f a p or the vector V a 
introduced below). 



3+1 geodesic equation and images in numerical spacetimes 



4 



Oe- To show it, we notice that the orthogonal projection P a := of p a on Oe's 

rest space is the linear 3-momentum of V with respect to Oe and ^ implies P l = EV 1 , 
which is exactly the relation between the 3-momentum and the 3-velocity of a particle 
(massive or not) of energy E. For a photon, the property p^ = 0, together with 
n M n M = —1 and fl8J), imply V^V^ = V^V 1 = 1, i.e. with respect to Oe, the photon travels 
at the speed of light (as it should!). For a massive particle, one has instead p^ < 0. 
The assumption E > 0, Q leads then to VjV % < 1: V cannot reach the speed of light. 

In the remaining part of this section, we do not consider a single geodesic, but a full 
congruence of them. This means that p a , E and V a are fields defined on the spacetime. 
Accordingly, we may consider their derivatives in any direction and not only along a 
geodesic as in Q. Let us then rewrite ^ by substituting ^ for p a and making use of 
(|2j and Q to express n^V p n a and V^n"; we get 

(n" + V^V^E (n a + V a ) + E (D a In iV + ra"V \y a - K^V + VV^V ) = 0. (9) 

Now, in the 3+1 formalism, the natural evolution operator is the Lie derivative £ m 
along the vector field m a := Nn a , for it preserves the property of being tangent to T, t 
[I]. Therefore we write 

n^V^V a = N^m^V^V* = N' 1 [£ m V a + V^V ^{Nn a )} 

= N- 1 £ m V a -K a /j y fl + V lM D lx \nNn a . (10) 
Similarly, since E is a scalar field, 

n^^E = N- l £ m E. (11) 
Also, the 4-dimensional and 3-dimensional covariant derivatives of V a are related by 

DpV a = rfsfpW, (12) 
from which we deduce the identity 

V^V„V a = V^D^V" - K^VV n a . (13) 

Note that, as in many places in the article, we have used the property n^V = along 
with expression Q for V/jn". 

Inserting (10), (11) and (13) into ([9]), we get, after division by E, 

N- 1 £ m V a + VD^V* - 2K a i y + E-\N- 1 £ m E + VD^V" + D a \nN 

+ [VDpkiN - K^VV + E-\N~ l £ m E + V^D fl E)] n a = 0. (14) 

Let us notice that the first line of this equation contains only terms tangent to £ t , 
whereas the term in the second line is manifestly parallel to n a . Hence the projections 
of ( 14 ) along n a and onto S t give respectively 

N- l £ m E + V^D^E + E{V^D^\nN -K^V V ) = (15) 
N~ l £ m V a + V^D^V* - 2K a /t V' i + E- 1 (N- 1 £ m E + VD^V" + D a \nN = 0. (16) 

These two equations involve only quantities intrinsic to Ef. We may then write them in 
a 3-dimensional form, using ^ to express the Lie derivative along m a : 

d 

£ m = 7^ - £p ■ (17) 
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If, in addition, we substitute (15) for £ m E in (16), we get 



\{jt~ £p ) E + VjDjE + EiyjD] lnN ~ K i kV3yk) = 

- £ p j V 1 + VWjV 1 - 27T j V j + V\K jk V j V k - V s D d In iV) 



N 

i /a 



+D i lnJV = 0. (19) 

This system constitutes the 3+1 geodesic equation in covariant form for a congruence 
of geodesies. 

2.3. 3+1 geodesic equation for a single geodesic 

Let us consider a specific member £ of the geodesic congruence, representing the 
worldline of a given particle V. In a coordinate system (x a ) = (t, x l ) adapted to the 



3+1 foliation (S t ) tgK (cf. Section 2.1), the equation of C can be written as 

x i = X\t), (20) 

where the X*'s are three smooth functions ffi. — > M.. This is nothing but a parametrization 
of C by the coordinate time t. Note that, a priori, t is not an afline parameter along C 
By definition, the velocity of V with respect to the Eulerian observer Oe is 

dt 

V% = ST' (21) 

dr E 

where dt is the displacement vector of Vs worldline with respect to the O^s one 
between t and t + dt, and dr E is the increment of Ce's proper time between t and t + dt. 
The origin of the (x l ) coordinates being "shifted" by the amount (3 l dt with respect to 
Ce's worldline, we have (cf. Figure [TJ 

dr E = iV dt and df = fidt + dX\ (22) 

Hence 

1 / \ r] 

V i = -[X i + ^ i ), with X l :=—. (23) 
The variation of E along C associated with the parametrization by t is 

where dj := d/dx j . Also, since E is a scalar field, £pE = f3 j djE and V j DjE = V j djE. 
) can be then written as 

E (NK jk V j V k - V j djN) . (25) 



Thanks to (|23j), (18) can be then written as 
dE 
~dt 

This evolution equation for the particle energy relative to the Eulerian observer is 
equivalent to equation (6) of Ref. [16J. The latter is actually an evolution equation in 
term of some affine parameter A along the geodesic, i.e. a parameter whose associated 
tangent vector is (up to some constant factor) the particle's 4- momentum: p a = dx a /dX. 
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Then, n^p^dX = n^dx^, which results in — EdX = —Ndt. Hence the relation between 
the two parametrizations of C: 

dX _ N 

dt ~ E' 



(26) 



Taking into account (26), we check that (25) is indeed equivalent to (6) of Ref. 

Let us now consider equation (19); we may express the Lie and covariant derivatives 
in terms of partial derivatives: 



9 £ 
N \dt~ £p 



V i + VWjV 1 



dV i 
N \~df 



.i'djV + V J <) r r ) + V'djV* + 3 r jk v s v 



where the r*- fc 's are the Christoffel symbols of the metric 7^ in E 4 . Then, by means of 
@ and the analog of (|24|) for V { : 

dV i 



dV l 



dt 



dt 



+ x 3 d*\r 



(27) 



we transform (19) into 



dV i 
~df 



NV j \V* (dj InN - K jk V k ) + 2K i j - 3 T) k V k ] - f j djN - V'djl?. 



Let us supplement this equation by the evolution equation for X 1 deduced from (23), 
to form the system 



dX 1 
~dT 

dV* 
~df 



NV i - ft 



(28a) 



NV j [V* (dj In N - K jk V k ) + 2K i j - 3 r jk V k ] - j^djN - \ V .(28b) 



Given the spacetime metric in 3+1 form, and hence the terms N, (3 l , j t3 \ 3 r* fc and 
Kij, (28) constitutes a system of six first order ordinary differential equations, that it 
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is sufficient to integrate with respect to t from initial data (X l (0),V l (0)) to get the 
geodesic worldline of V in the form (20). 



Note that, contrary to ( 18 )-( 19 ), the energy equation (25) and the system (28) 
are meaningful for a single geodesic: they involve only the derivatives dE/dt, dX l /dt 
and dV l /dt, which are derivatives along the geodesic (and not transverse to it). To 
strengthen this point, we present in Appendix A an alternative derivation of ( p5| and 
(28), which does not rely on the assumption that C belongs to some geodesic congruence. 

A 3+1 form of the geodesic equation has already been derived by Hughes et al. [32] 



(cf. also Section 7.2 of the textbook [9]). However, it differs from the present one by 
the following features: (i) it involves the components Pi of the 4-momentum instead of 
V 1 , (ii) the evolution parameter is the affine parameter A and not the coordinate time t, 
(iii) it is valid only for massless particles. Moreover, in Ref. [TJ] no evolution equation 
for the particle's energy, equivalent to our equation (25) is provided. Note also that the 
3+1 geodesic system of [12] is applied to the determination of the event horizon, not to 
the formation of images. 

In |Appendix B we combine (28a) and (28b) into a single second order equation for 
X l (t). We recover in this way the standard 4-dimensional geodesic equation. However, 
this equation is less convenient for numerical integration for it involves time derivatives 
of the lapse and shift, contrary to the system (28). 



3. Redshift factor 

3. 1 . General formula 



The integration of (25 ) forward in t gives the energy of the particle V with respect to the 
Eulerian observer at any point. Let us consider the case in which V is a photon emitted 
at some event A by an observer O em (the "emitter") and received at some event B by 
an observer O rec (the "receiver" ) . Note that these observers are not necessarily Eulerian 
observers. In this way, the problem considered here generalizes that of Ref. |16| . which 
was limited to Eulerian observers. The redshift factor z is defined by 

l + z = ^, (29) 

^rec 

where v era (resp. v rec ) is the photon frequency measured by O em (resp. O mc ). The 
frequency being related to the energy by the Planck-Einstein formula e = hu, the above 
relation can be written 



e 



where e em = — P^Wem (resp. e rec = — Py\ B u^: ec ) is the photon energy with respect to 
C em (resp. O rcc ) and u® m (resp. u" ec ) is the 4- velocity of O cm at A (resp. of O rec at B). 
Let us perform the 3+1 decomposition of these 4- velocities: 

< m = r em (n a + U:j with n„EC = (31) 
< ec = r rec (n Q + £C) with VC = 0. (32) 
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r em (resp. r rec ) is then the Lorentz factor of O em (resp. O rcc ) with respect to the 
Eulerian observer C E and U^ m (resp. U® c ) the 3-velocity of O em (resp. O rcc ) with 



respect to Oe- From the normalization relation 



(1 - jijU^U^) 



-1/2 



and T r 



Combining (\8v and (31) yields 



E\ A T em (l- ltJ :V' 



— 1, we get 

(i-^uLuL) 



u j ) 

em/ 



-1/2 



(33) 



(34) 



A similar relation holds for e rec , so that (30) becomes, once (|33j) is taken into account, 

1/2 

i + : = ^r 1 - — i 1 — "-';, lv, ;:;" : I . (35) 



E\ A l- lljV *\ A U im 

E\ B i- lijV t\ B uL 



The procedure to compute the redshift factor is then as follows 
date[[| E \a = E ( t = Ia), V*\ a = V\t = t A ) and X\t = t A ) 



Given some initial 



x 



A- 



where 



.r 



are the 



coordinates of the event A on the hypersurface one integrates the system formed 
by pHb and (pb to get E\ B = E{t = t B ) and V'% = V\t = t B ). Then, one uses 



(35) to evaluate z. Note that since (25) is a homogeneous equation in E and only the 
ratio E\ A j E\ B is involved in the expression of z, the initial value E\ A can be chosen 
arbitrarily. Of course, if one wants to manipulate some physically relevant value of E, 
one may deduce E\ A from the photon energy with respect to the emitter, e era , via (34). 



3.2. Limiting cases 



As a check of formula (35), let us consider the special case of an inertial observer in 



Minkowski spacetime receiving a photon emitted by a moving source. Choosing for 
(Ej)tgK the time foliation associated with that inertial observer, we have Oe = O 
so that Ul ec = 
implying E\ B 



0. Moreover, N = 1 and = 0, so that (25) reduces to dE/dt 



rec; 

0, 



E\ A . Accordingly (35) reduces to 



l + z 



1-/, 



V 1 1 TP 

ij v I A em 



1 fij ^em 



UL 



(36) 



where fij is the flat metric. We recover the special relativistic formula for the Doppler 
effect. In particular, if the photon travels in the same direction (up to a sign) as the 

U 2 with U : 



emitter, we have fij V i \ A U£ m = U and fijU^JJ^ 



± 



fijU* m Ul m , with 



a + (resp. — ) sign if the emitter is approaching to (resp. receding from) the receiver, 
so that ([361 gives the well-known formula 



l + z 



1 - U 
l + U' 



(37) 



Another check of formula (35) and equation (25) is provided by the propagation 
of a photon between two static observers in Schwarzschild spacetime. Choosing (E t ) 

|| In a ray-tracing code, the integration is usually performed backward, i.e. from B to A. Accordingly, 
the roles of A and B have to be swapped in the following discussion. 
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to be the standard foliation associated with Schwarzschild time coordinate t, the static 
observers coincide with Eulerian observers, so that we have U* m = and U* ec = 0. 



Accordingly, (35) reduces to 



I : = |}a (38) 



On the other side, (25) reduces to 



f = (39) 

because = for the foliation (£ t ) and (5 l = in standard Schwarzschild coordinates 
(t, r, 6, if). Since dN/dt = 0, we may rewrite the above equation as 
dE EdN 

M = -NW (40) 
from which we deduce immediately 

±(EN) = 0. (41) 



Hence E = const/iV and (38) becomes 1 + z — N\ B / N\ A . Using the value of the 
lapse in terms of the mass parameter M of Schwarzschild metric and the Schwarzschild 
coordinate r, N = yl — 2M/r, we get 



We recognize the classical formula for the gravitational redshift (Einstein effect). In 
particular, in the "Sirius B configuration" (re — > +oo and M/ta «C 1), we get 
z ~ M/ta, as it should be. 

4. Applications 

4-1. Implementation in the GYOTO code 



The 3+1 geodesic equations (25 ) and (28 ) have been implemented in the ray-tracing code 
GYOTO pm ITTj . The integration in t is performed by means of a fourth-order Runge- 
Kutta algorithm. The 3+1 fields (A, (3\ 7y, Kij) have to be provided by an external 
code. An example of using the 3+1 fields from a numerical relativity spectral code is 
given in Section 3 of [TO], where it is shown how to get the values of (N, (3 l ,7^, K^) at 
each point of the geodesic from the outputs of the spectral code. 

The derivatives with respect to some affme parameter A along the 
geodesic, (dt/dA, dr/dX, d9/d\, d(p/dX), can be derived from the 3+1 derivatives 
(dr/dt, dd/dt, dtp/dt) provided that one knows the value of dt/dX. The latter is de- 



duced from the value of E resulting from the integration of (25), by noticing that 
E = —n^p^ = N p l and p t = mu l = mdt/dX for a massive particle (A is then the parti- 
cle's proper time) and p l = dt/dX (up to some constant rescaling of A) for a photon. 
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Figure 2. Relative error, in units of 10 -5 , on a null geodesic integrated by GY0T0 in 
a numerically computed Kerr metric with spin parameter a = 0.5 M, when compared 
to the standard integration using the analytical expression of the Kerr metric, for 
the Boycr-Lindquist coordinates r (yellow), 9 (green) and ip (magenta) as a function 
of time coordinate t. The geodesic is integrated backward in coordinate time from 
t = 1000 M,r = 100 M until t = 0,r = 865 M. It reaches the smallest distance 
(r = 4.3 M) from the black hole around t — 900 M: this is where the error is the 
largest. 



4-2. Numerical tests 

A preliminary test of the 3+1 computation of geodesies in a numerical spacetime has 
been provided in Figure 8 of pU], where a timelike geodesic computed by the 3+1 
method was compared to that computed by integrating the standard 4-dimensional 



geodesic equation [equation (B.l) below]. The spacetime was that of a rapidly rotating 
relativistic star, numerically generated by means of the LORENE/nrotstar code [TTJUH]- 
We present here a more detailed test, regarding a null geodesic around a Kerr 
black hole, with a spin parameter a = 0.5 M (M being the black hole mass). The Kerr 
spacetime is described in Boyer-Lindquist coordinates (i, r, 9, <p) and its 3+1 "numerical" 
version has been prepared on a spectral grid via a code using the L0RENE library [18J. 
A test null geodesic, that comes close to the event horizon and therefore subject to 
strong-field effects, has been integrated by GY0T0 via two methods: (i) integration 
of the 3+1 geodesic equations in the numerically generated Kerr spacetime, and (ii) 
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standard 4-dimensional integration using the analytical expression of the Kerr metric in 
Boyer-Lindquist coordinates [TO] . 

Figure [2] shows the resulting relative difference between the spatial coordinates 
(r(t),0(t),<p(t)) obtained by the two methods. The maximum relative error, occurring 
when the geodesic comes at the closest distance to the black hole, is of a few 10 -5 , which 
is very satisfactory. 

4-3. Images of a stationary rotating neutron star 

Using the 3+1 geodesic equations implemented in GYOTO, we have computed the image 
perceived by a distant observer of a rapidly rotating neutron star in a spacetime 
computed by the LORENE/nrotstar numerical code [T7J [J_2]. The neutron star model is 
built upon an equation of state derived by Akmal, Pandharipande and Ravenhall [19] 
and is described in detail in Section 3.5.3 of [T7]. The mass of the star in M = 1.4M 
and it is chosen to be either non-rotating or rotating at the frequency of 716 Hz (the 
largest observed frequency |20j). 

Figure [3] shows the image of these two models of neutron stars, assuming the surface 
of the star is optically thick and emits as a black body at 10 6 K. The effect of relativistic 
beaming, due to the star's rotation, appears clearly on the right panel. Moreover, the 
rotating star is oblate, the ratio of its apparent polar radius to its apparent equatorial 
radius is 94%. This effect is a known consequence of its rotation. 

Since the spacetime is stationary and axisymmetric, two quantities must be 
conserved along each geodesic: the components pt and p v of the 4-momentum. A third 
constant of motion is the squared norm of the 4-momentum : p^ = (null geodesic). 
The constancy of these three quantities is not imposed in the code. We therefore monitor 
them along the geodesies in order to check that the integration is performed correctly. 
In the present case, the squared norm of the photon p^ stays below a few 10~ 5 , the 
maximum relative error on p t is a few 10~ 6 and the maximum relative error on p v is a 

few irr 4 . 

Let us end this section by a remark that concerns also the next section. Considering 
a neutron star at a distance of 1 kpc with a size of 10 km, its apparent size disregarding 
any relativistic effect on the photon's trajectory would be less than 10 -10 arcsecond. This 
is of course far beyond the resolution of any current or near future instrument. The 
present and next sections must therefore not be read as describing possible observational 
tests. 

4-4- Images of the collapse of a neutron star to a black hole 

To illustrate the 3+1 geodesic computation in dynamical spacetimes, we consider the 
astrophysical scenario of an unstable non-rotating neutron star collapsing to a black 
hole. This scenario is numerically modeled using the CoCoNuT code [21] , which solves 
the relativistic hydrodynamics equations, coupled to the Einstein equations for the 
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Figure 3. Images (i.e. map of specific intensity) of a non-rotating (left) and 716 Hz 
rotating (right) stationary neutron star, with an optically thick surface emitting black 
body radiation at 10 6 K. The color bar is common to the two panels and is given in 
SI units, Wm~ 2 ster -1 Hz -1 . The frequency of the photons in the observer's frame is 
chosen to be 10 17 Hz, close to the maximum of the Planck function at 10 6 K. 



gravitational field, within the so-called conformal flatness condition (CFC). In the multi- 
dimensional case, CFC is an approximation to general relativity where the 3-metric 



where ip is the conformal factor and a flat 3-metric. In our particular case of 
spherical symmetry (no rotation), this is not an approximation but reduces to the choice 
of isotropic coordinates. Therefore, the whole non-rotating simulation can be exactly 
described within CFC, even after the formation of the black hole's apparent horizon. 

Initial models of spherical neutron stars are equilibrium configurations on the 
unstable branch, computed in isotropic gauge using the LORENE/rotstar_dirac 
code [TBI E2]. The equation of state used for generating initial data and for computing 
the collapse is a polytropic one, neglecting any temperature effects; pressure p is related 
to baryon density n through 



with the adiabatic index 7 = 2 and k = 4.01 x 10~ 56 in SI units. 

The neutron star used as initial data for the CoCoNuT code is a 1.62M (gravitational 
mass) configuration, corresponding to a 1.77 Mq baryon mass. Its central density is 
1.56 fm~ 3 , central lapse N c = 0.40 and circumferential equatorial radius -R c j rc = 10.4 
km. The global accuracy indicator gives 4 x 10~ 9 . To this configuration, we add a 
perturbation to the density profile ensuring that the unstable star collapses to a black 
hole and does not evolve to the stable branch (migration). The density profile is modified 
according to 



7i 3 - (pj) is conformally flat: 

7y ^ fiji 



(43) 



p = ku 



7 



(44) 




(45) 
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where r is the coordinate radius, A = 0.01 is the relative amplitude of the perturbation 
and Ro = 10 km its typical size. For the CoCoNuT code, we use 500 radial cells on a 
uniform grid. 

The collapse proceeds as expected (see e.g. [23J), until the formation of an apparent 
horizon detected by the finder described in [24], at t = 0.438 ms after the beginning 
of the collapse. The simulation is stopped at t — 0.495 ms, when all the matter has 
entered the black hole (up to numerical accuracy). The run is stopped because of the 
too strong increase of the gradients in many quantities (e.g. the conformal factor ip) 
near the center of the star. This is due to the use of maximal slicing gauge condition 
(trace of set to zero), which has the well-known property of yielding a singularity- 
avoiding time slicing. Nevertheless, as stated before, most of the matter has entered 
the black hole at that time and there is no longer evolution of the black hole. During 
the computation of the collapse, quantities which are used to integrate the system (28) 
are exported from CoCoNuT to GY0T0 following the procedure described in Section 3.3 
of [10]. These quantities are the 3+1 metric and related fields (N, fl 1 ,"/^, Ky), together 
with the fluid 4- velocity u^ uid , the radius of the neutron star and the location of the 
black hole apparent horizon. 

When integrating a null geodesic, GY0T0 uses an interpolation at third order in 
the time coordinate to determine the value of the 3+1 fields at each integration step. 
Each geodesic is integrated backward in time until it either reaches the star's surface, 
or the black hole's event horizon. The difficulty here is that the location of the event 
horizon is not known by CoCoNuT, only that of the apparent horizon (that lies inside the 
event horizon) is known. The integration of geodesies that reach the star after the event 
horizon radius has become larger than the star's radius is thus non trivial. There are 
thus two stop conditions for geodesies that reach the central object (i.e. for geodesies 
not escaping towards infinity): 

• the star's surface is hit. In such cases, the specific intensity emitted by the hit point 
can be computed. 

• the fourth order Runge-Kutta adaptive step becomes smaller than a given lower 
limit (fixed to 10~ 6 M). This latter case corresponds to a geodesic "accumulating" 
near the event horizon. Let us remind the reader that a backward integrated 
geodesic can never cross the event horizon (by the very definition of an event 
horizon) . 

The value of the integration step lower limit is chosen in such a way that the image of the 
event horizon on the observer's screen is smooth and spherically symmetric (choosing a 
too big limit results in a non spherically symmetric, noisy event horizon image). 

Moreover, the norm of the 4-momentum is not conserved in the very last integration 
steps close to the apparent horizon. This is due to the fact that the quantity dt/dX 
becomes very large (it should, actually, diverge if the integration were perfect) at the 
event horizon: as the norm is proportional to this quantity, it is no longer conserved. 
As a check of this fact, the evolution of the norm divided by dt/d\ has been considered. 
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Figure 4. Images (i.e. map of specific intensity) of a non-rotating collapsing neutron 
star, with an optically thick surface emitting black body radiation at 10 6 K. The color 
bar is common to the four panels and is given in SI units, Wm -2 ster -1 Hz -1 . The 
frequency of the photons in the observer's frame is chosen to be 10 17 Hz, close to the 
maximum of the Planck function at 10 6 K. 

This new quantity stays close to zero even at the horizon. 

Figure [4] shows four images of a collapsing non-rotating neutron star, as perceived 
by a distant observer. The first image is computed before the start of the collapse: it 
is thus the image of a stationary (unstable) neutron star. The other three images show 
different stages of the collapse, until the whole star nearly disappears below the event 
horizon. The intensity is shown in logarithmic scale, since the very high gravitational 
redshift leads to a very high dynamic range in the last images. The event horizon first 
appears at the center of the star due to the fact that this part of the star is closer to 
the observer: photons at different parts of the image have been emitted at different 
coordinate times t, later times for the central parts of the image, earlier times for the 
external parts of the image. As a consequence of this fact, the coordinate radius of the 
star at which a given photon is emitted is not the same for all pixels in a given image: it 
is shorter at the center of the image (more evolved part) than at the edge. For instance, 
the coordinate radius of the star at the emission of the photon reaching the central pixel 
of the second image on the left is of 2.7 km while it is 4.7 km on a pixel located at the 
edge of the star's image. The coordinate radius of the star on the left image is of about 
7 km (here, the coordinate radius is the same at any pixel on the image, as the left 
star is stationary, and smaller than the circumferential radius cited above of 10.4 km), 
whereas the coordinate radius of the star in the rightmost image is of approximately 
2.9 km for a pixel at the edge of the image. The fact that the ratio of the apparent radii 
of the star on the left and right panels of Figure [4] is much less than 7/2.9 ~ 2.4 is due 
to the strong bending of geodesies in the vicinity of the nascent black hole, resulting 
in its apparent radius being larger as the object becomes more compact. This effect is 
exactly the same as the one which makes the black hole shadow (the black area in the 
image of a black hole in front of an emitting region) larger than the projected size of 
the event horizon. For a Schwarzschild black hole, the enlarging factor is 3^f3/2 ~ 2.6. 
This explains why the shrink of the size of the image of the collapsing star shown in 
Figure [| is not very pronounced. Finally, let us stress the fact that the time elapsed 
between the appearance of the event horizon and the disappearance of the whole star 
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behind it is extremely short. It is of approximately 0.2 ms in the observer's frame. 

Figures [3] and [4] are the first examples of ray-traced images in numerically computed 
spacetimes. In the present work, the physics of emission at the surface of the neutron 
star has not been studied in detail at all. In particular, the emission during the collapse 
will most certainly be quite different from a simple blackbody. However the aim of the 
computations presented here is not (yet) to propose astrophysically relevant images, 
but to give first examples of the interest of such a ray-tracing code as GY0T0, capable of 
integrating geodesies in numerical spacetimes. Future works will be devoted to applying 
this code to diverse astrophysically realistic scenarios. 



5. Conclusion 



We have re-expressed the geodesic equation within the framework of the 3+1 formalism 
of general relativity, obtaining equations (25), (28) and (35). Equation (25), ruling 
the evolution of the particle energy with respect to Eulerian observers, has already 
been derived (in an equivalent form) by Merlin and Salgado [16]. On the other hand, 
the system ( [28] ) for the position of the geodesic and the redshift formula (35) are 
novel. In particular, (28) significantly differs from previous 3+1 geodesic equations 
in the literature [12], as discussed in Section 2.3[ The 3+1 equations have been 
implemented in the ray-tracing code GY0T0 [TU1 [TTJ, which has enabled us to compute 
images of stationary and collapsing neutron star numerical spacetimes generated by the 
LORENE/nrotstar [T7J [18] and CoCoNuT [21] codes. 

Future work will be devoted to the development of ray-tracing computations in 
numerical spacetimes for astrophysically relevant problems. In particular, we shall try 
to carry on computations of images and spectra of astrophysical phenomena in the 
vicinity of compact objects, which can be alternative to black holes. This work is of 
particular interest in the perspective of a direct test of the nature of the central compact 
object of the Galaxy, Sgr A* (see the review [26j, and in particular the section devoted 
to the alternatives to the black hole case). 

The capability of GY0T0 to integrate geodesies in numerical spacetimes will be very 
interesting too, in order to visualize spacetimes, be it binary black holes spacetimes, 
binary neutron stars spacetimes, black hole - neutron star binary spacetimes, or any 
other interesting metric (see reviews on these topics by [271 EHJ EHJ [30]) . GY0T0 could 
be used to image a background sky of stars, or a simple coordinate grid, putting in light 
the effect of strong gravity on these background objects. 
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Appendix A. Derivation for a single geodesic 



Here we consider a single (null or timelike) geodesic £ and not a congruence as in 



Sect. |2.2| In the context of the 3+1 formalism, a natural parameter along £ is the time 
coordinate t that labels the foliation (£ t ) te R. A priori, t is not an affine parameter along 



£; it is related to the affine parameter associated to p a by (26). The vector tangent to 
£ associated to the parametrization by t is q a := dx a /dt, dx a being the infinitesimal 
displacement vector along £ corresponding to the infinitesimal parameter increment dt. 
The vector q a obeys q^V^t = 1. Moreover, q a and p a being vectors both tangent to £, 
they must be collinear. We deduce then from pi) that 



Q 



N(n a + V a ). 



(a.i; 



Here, the fact that the proportionality coefficient in (A.I) is N comes from g M V at 



n^V^t = N' 1 and V*V„t 



0. Taking into account l\8h, the geodesic equation fl7h is 



equivalent to 



g^V M [E{n a + V a )\ = 0. 



(A.2) 



Expanding and using Q, as well as (A.I), we get 

q^V^E (n a + V a ) + E (D a N - NK^V + q^V^V ) = 0. (A.3) 

Note that the only derivatives of E and V a that appear in this equation are derivatives 
along £ (through the operator q^V M ). Consequently (A.3) is valid for a single geodesic, 
contrary to (fl4l) which holds only for a congruence. 



The projection of (A.3) along n a yields 
-q^W^E + Eq^V^V" = 0, 
where we have used n u V u = 0, n u D u In N = and n v K v 



0. Now, since n,,V L/ = 0, 



we have n u V \y v = -V U V M n„, i.e., thanks to Q, n v V ' = V v {K tlv + D u InNn^). In 
addition, from the very definition of the vector q a , q^V^E = dE/dt. We thus end with 

dE 



dt 



EV j (NK jk V k - djN) 



(A.4) 



which is exactly (25) 



A" 



The orthogonal projection of (A.3) onto E t is performed via the operator 7 
Thanks to the properties "i a v n v = 0, ~f a v V v = V a , ^ a v D v N = D a N and Tu K % = K 
this yields 

<f V M £ V a + E (D a N - NK a ^V>* + i a v q^V ^V u ) = 0. (A.5) 



To evaluate the last term, let us expand the velocity vector V a onto the coordinate basis 
ifia) '■= (d/dx a ) associated to the coordinates (x a ) = (t,x l ). Since V a is tangent to S t 
it has no component along eo = d/dt: 
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We then have, since •y a u e^ = (for e 3 - is tangent to E t ), 

7 a /V,V = <fV^ e« + V'-y^Vrf 

= g"V M V J ' e« + WV„(n" + ^)V M e} 

= g"V^' e* + W j (fXV.eJ + VD^efj . 
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(A.7) 



Now, thanks to (10) with V a replaced by e", we can write 

7"/V^ = 7% (A^™ - + In iV n") 



= N- 1 £ ei p a -K a ll e$, 
where we have used (17), r y a u K u ^ = K a „, r y a v n u - 



(A.8) 



[eo,ej] = (since eo and 
ej are vectors of a coordinate basis), £pe°- = —£ ej (3 a , and ■y a u £ ej (3 U = £ ej (3 a (since 



0? £e e j 



e" and (3 a are both tangent to S t ). Thanks to (A.8), (A.7) becomes 

7VV^" = q"V»V j e* + ^ [£ e . a + N (VDrf - K a ^)] . (A.9) 

Substituting this expression into (A. 5) and setting a — i (for a = 0, the equation 
reduces to = 0), we get 

<f V M l/ j e) = -E-^VpE V l V j [N (K* k (Jj - V k D k e)) - £ e . ft] - D l N + NK^VK 
Let us consider the components of this vector equation in the coordinate basis (ej). 



We then have e* 



D k e) 



3r« 



jki 



D l N = Y^djN and, by the definition of a Lie 
derivative, £ ej ft = djft. Since in addition, q^W tl V j = dV { /dt and q^V^E = dE/dt, 
we get 

— = -^-^ V i + NV j (2K l j - 3 T l jk V k ) - ft j d 3 N - V j djft. (A. 10) 
hi dt " 



dt 



Substituting (A. 4) for dE/dt, we recover (28b). 



Appendix B. Second order form of the 3+1 geodesic equation 

The standard (4-dimensional) geodesic equation is 
d 2 X a +4pa cU^cLY" _ Q 



dA 2 



dA dA 



(B.r 



where (i) X := t, in addition to the X"s defined by (20), (ii) the T a 's are the 



Christoffel symbols of the spacetime metric g a p and (iii) A is an affine parameter along 
the particle's worldline C If £ is timelike, A is equal to a constant times the particle's 
proper time r. More specifically, if A is the affine parameter associated to the particle's 
4-momentum and m is the particle's mass, we deduce from the relations p a = dx"/dA 
and p a = mdx a /dr that A = r/m. 



Let us check that (B.l) can be recovered from the system (28). Extracting V 



from (28a), substituting it in (28b) and expanding, we obtain a second order differential 
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equation for X l (t): 
1 



X l + 



N 



K jk (X j + (3 3 )(X k + (3 k ) - (2X3 + pj) djN 



Dj? - NK l 3 + ^(K jk (3 k - djN) 



+Nj %3 djN - 2NK l j f3* + '- [ K jk p>p 



dN 
~dt 

& + 

dN 
~ ~dt 



X 1 



3 r , t + ?-K ik ) Px k 



■>'•- X 

P 3 djN 



dl3 i 

+-^r + P'Dj? = 0. 



Of 



(B.2) 



On the other hand, the 4 r" 's appearing in (B.l) can be expressed in terms of the 



3+1 quantities as follows (cf. Appendix B of [8]): 
1 (dN 



4-pO 



00 



4-pO 



0j 



+ PfyN - K jk p j p k 
(d 3 N - K jk p k ) 



N V dt 
1 

N 



4-pO 



N^ 



:K 



(B.3) 

(B.4) 
(B.5) 



4 ni 



P* 



00 



ON 



N^djN - 2NK 1 ^ + tL ( K jk p>p k - — - ptfyN ) + + pWjP* (B.6) 



dt 



dp 1 



Of 



4-pi 



P* 



0j 



4pi 



D i P i -NK i j + ^(K jh p k -d j N) 



In addition, we have 



dX» _ . dt 
dA dA 



d 2 X a 



x c 



dt 

dA 



+ X< 



eft 

dA 2 ' 



(B.7) 
(B.8) 

(B.9) 
(B.10) 



Accordingly, for a = 0, (B.l) becomes (note that X° = 1 and X° = 0): 
d 2 t 



dA 2 



+ 



dt 
dA 



4 r° 00 + 2 4 r° 0j x^ + ^ jk Px k 



In view of (B.3)-(B.5), we get 





- 2 d 2 t 


1 


(9 


dA 2 





K jh (X> + P j ){X k + P k ) - (2X3 + pj) djN _ 



0. 



dN 
~dt 



(B.ll) 



For ct = i, pi] ) becomes, thanks to (|B79|)-(|B.10|) 

dt 



X l + 



"2 J2 



d 2 t - 



dA ) dA 2 



a* + 4 r 00 + 2 4 r n , a j + 4 r „ a j A fc 



0.r 



In view of (B.ll) and (B.6)-(B.8), we recover (B.2). 
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